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1 .  INTRODUCTION 


The  finite  element  analysis  of  general  shell  structures  has  been 
a  very  active  field  of  research  for  a  large  number  of  years  [14,  29]. 
However,  despite  the  fact  that  many  different  shell  elements  have  already 
been  proposed,  the  search  for  a  shell  element  capable  of  representing  the 
general  nonlinear  behavior  of  shells  with  arbitrary  geometry  and  load¬ 
ing  conditions  in  an  effective  and  reliable  manner  is  still  continu¬ 
ing  very  actively. 

During  recent  years  it  has  become  apparent  that  two  approaches  for 
the  development  of  shell  elements  are  very  appropriate: 

-  The  use  of  simple  elements,  based  on  the  discrete-Kirchhoff 
approach  for  the  analysis  of  thin  shells  [2,  5  -  9]. 

-  The  use  of  degenerated  isoparametric  elements  in  which  fully 
three-dimensional  stress  and  strain  conditions  are  degenerated 
to  shell  behavior  [2,  3,  5,  7,  17,  19,  24,  29], 

The  latter  approach  has  the  advantage  of  being  independent  of  any 
particular  shell  theory,  and  this  approach  was  used  in  ref.  [3]  to 
formulate  a  general  shell  element  for  geometric  and  material  nonlinear 
analysis.  This  element  has  been  employed  very  successfully  when  used 
with  9  or,  in  particular,  16  nodes.  However,  the  16-node  element  is 
quite  expensive,  and  although  it  is  possible  to  use  in  some  analyses 
only  a  few  elements  to  represent  the  total  structure  (see  Sections  4.2 
and  4.5),  in  other  analyses  still  a  fairly  large  number  of  elements 
need  be  employed  [5], 

Considering  general  shell  analyses,  much  emphasis  has  been  placed 
onto  the  development  of  a  versatile,  reliable  and  cost-effective 
4-node  shell  element  [16,  17,  22,  28].  Such  element  would  complement 
the  above  high-order  16-node  element  and  may  be  more  effective  in 
certain  analyses.  The  difficulties  in  the  development  of  such  element 
lie  in  that  the  element  should  be  applicable  in  a  reliable  manner  to 
thin  and  thick  shells  of  arbitrary  geometries  for  general  nonlinear 
analysis . 

The  objective  in  this  paper  is  to  present  a  simple  4-node  general 
shell  element  with  the  following  properties: 

-  the  element  is  formulated  using  three-dimensional  stress  and 
strain  conditions  without  use  of  a  shell  theory; 

-  the  element  is  applicable  to  thin  and  thick  shells  and  can  be 
employed  to  model  arbitrary  geometries; 
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-  the  element  is  applicable  to  the  conditions  of  large  displace¬ 
ments  and  rotations  but  small  strains,  and  can  be  used  effec¬ 
tively  in  materially  nonlinear  analysis. 

The  formulation  of  the  element  is  quite  simple  and  transparent, 
and  the  element  has  good  oredictive  capability  without  containing 
spurious  zero  energy  modes. 

In  the  next  section  of  the  paper  we  discuss  some  basic  considera¬ 
tions  with  respect  to  the  assumptions  used,  and  in  Section  3  we  present 
the  element  formulation  for  nonlinear  analysis.  The  results  obtained 
in  numerical  solutions  that  demonstrate  the  properties  of  the  element 
are  given  in  Section  4. 


2.  BASIC  CONSIDERATIONS 


The  formulation  of  the  4-node  shell  element  represents  an 
extension  of  the  shell  element  discussed  in  refs.  [2,  3],  and  we 
therefore  use  the  same  notation  as  in  those  references.  Also,  to  focus 
attention  onto  some  key  issues  of  the  formulation,  we  consider  in 
this  section  only  linear  analysis  conditions. 


The  geometry  of  the  element,  see  Fig.  1,  is  described  using 
[2,  p.  255]: 


4 


4 


(1) 


k=l 


k=l 


where  the  are  the  two-dimensional  interpolation  functions 

corresponding  to  node  k;  the  r.  are  the  natural  coordinates;  and 


0 

x.  =  Cartesian  coordinates  of  any  point  in  the  element; 

o  k 

x..  =  Cartesian  coordinates  of  nodal  point  k; 

0  k 

V  .  =  Components  of  director  vector  at  node  k  (which  is  not 


necessarily  normal  to  the  midsurface  of  the  element); 


B 


r  2 


Figure  1  Four-node  shell  element 
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•  I ?  k 

and  a,  is  the  shell  thickness  at  node  k,  measured  along  the  vector  V 


The  left  superscript  is  zero  for  the  initial  qeometry  of  the  element  and 
is  equal  to  1  for  the  deformed  element  geometry.  Note  that  the  thickness 
of  the  element  varies  and  the  element  is  in  general  non-flat. 

The  displacements  of  any  particle  with  natural  coordinates  r.  of  the 
shell  element  in  the  stationary  Cartesian  coordinate  system  are: 


4 


4 


k=l 


k=l 


L' 

where  the  u.  are  the  nodal  point  displacements  into  the  Cartesian 
coordinate  directions,  and  the  and  6^  are  the  rotations  of  the 


director  vector  k  about  the  ^V,k  and  axes  (see  Fig,  1). 
— n  —  i 


A  basic  problem  inherent  in  the  use  of  the  above  interpolation 
of  the  displacements,  and  the  derivation  of  the  strain-displacement 
matrices  therefrom,  is  that  the  element  "locks"  when  it  is  thin.  This 
is  due  to  the  fact  that  with  these  interpolations  the  transverse  shear 
strains  cannot  vanish  at  all  points  in  the  element,  when  it  is  subjected 
to  a  constant  bending  moment.  Hence,  although  the  basic  continuum 
mechanics  assumptions  contain  the  Ki rchhoff  shell  assumptions,  the 
finite  element  discretization  is  not  able  to  represent  these  assumptions, 
rendering  the  element  not  aDplicable  to  the  analysis  of  thin  plates  or 
shells  (see  [2,  d.  240]  and  [5,  7]).  To  solve  this  deficiency,  various 
remedies  based  on  selective  and  reduced  integration  have  been  proposed 
[1  7,  22,  23]  but  there  is  still  much  room  for  more  effective  and  reliable 
elements  for  general  nonlinear  analysis. 

Considering  our  element  formulation  —  because  the  problem  lies  in 
the  representation  of  the  transverse  shear  strains  —  we  proceed  to  not 
evaluate  these  shear  strains  from  the  displacements  in  Eq.  (2)  but 
to  introduce  separate  interpolations  for  these  strain  components. 

Since  we  consider  non-flat  shell  elements,  the  separate  i nt&*molations 
are  performed  effectively  in  a  convected  coordinate  system.  ' 


^Note  that  in  refs.  [2]  and  [3],  the  shell  element  formulation  is  dis¬ 
cussed  in  the  global  stationary  coordinate  system,  because  all  displace¬ 
ment  components  are  interpolated  in  the  same  way.  To  emphasize  that 
we  use  here  stress  and  strain  measures  in  the  convected  coordinate 
system,  we  place  a  curl  (~)  over  these  quantities. 
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The  choice  of  the  interpolation  for  the  transverse  shear  strain 
components  is  the  key  assumption  in  our  element  formulation,  because 
adequate  coupling  between  the  element  displacements  and  rotations  must 
be  introduced  and  the  element  should  not  exhibit  any  spurious  zero 
energy  modes.  For  our  element  we  use,  see  Fig.  2, 


e13  =  \  (1  +  r2}  £13  +  I  d  ~  r2}  gi3 
e23  =  (i  +  rx)  i^3  +  \  d  -  ri)  £33 


(3) 


Since  the  kinematic  relations  for  the  above  shear  strains  are  not 
satisfied  using  Eq.  (3),  we  impose  them  using  Lagrange  multipliers 
[2,  27]  to  obtai n 


n* 


)  dV 


V 


+ 


,23  r 
X  (e 


23 


~DI  \  ... 
e23>  dV 


-  w 


(4) 


where  the  t1J  are  the  contravariant  components  of  the  Cauchy  stress 

tensor  [13,  15],  the  £ ■  ■  are  the  covariant  components  of  the  infini- 

'J  13  ?3 

tesimal  strain  tensor,  the  X  and  X  are  the  Lagrange  multipliers, 

the  and  £^3  are  the  transverse  shear  strains  evaluated  using  the 

displacement  interpolations  in  Eq.  (2),  and  W  is  the  potential  of  the 

external  loads.  For  the  Lagrange  multipliers  we  choose  the  following 
interpolations , 


11 


.r3 


e^3  interpolation 


Figure  2  Interpolation  functions  for  the  transverse 
shear  strains 
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(5) 


X13  =  XA  fi(rj)  6(1  -  r2)  +  XC  6(r1)  6(1  +  r2) 


X23  «  XD  6(r2)  6(1  -  rx)  +  XB  6(r2)  6(1  +  r}) 

where  6(...)  is  the  Dirac-delta  function.  This  represents  a  weakening 
of  the  Lagrange  multiplier  constraint  in  Eq.  (4)  [10].  Substituting 
from  Eq.  (5)  into  Eq.  (4)  and  invoking  that  6n*  =  0  gives  the  distinct 
constraints 


.01 

e13 


at  A 


at  D 


.01 

£13 


at  C 


(6) 


at  B 


Hence,  the  complete  element  stiffness  matrix  is  calculated  using  the 
functional 


(7) 


with  stress  and  strain  components  in  convected  coordinates  and 

-  Eqs.  (1)  and  (2)  to  evaluate  the  strain  components  e22 
and  e12; 

-  Eq.  (3)  to  evaluate  the  strain  components  i^y  and 

-  Eq.  (6)  to  express  the  variables  ^3*  *^3  ancl  ^23  in 
terms  of  the  nodal  point  displacements  and  rotations  of  Eq.  (2). 
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Considering  the  reDresentation  that  we  have  chosen  for  the  trans¬ 
verse  shear  strains,  we  can  make  the  following  three  important  obser¬ 
vations  : 

1)  The  element  is  able  to  represent  the  six  rigid  body  modes. 

The  element  contains  the  rigid  body  modes  because  zero  strains  are 
calculated  in  the  formulation  when  the  element  nodal  point  displacements 
and  rotations  correspond  to  an  element  rigid  body  displacement.  This  can 
be  verified  by  using  Eqs.  (1)  to  (6)  to  evaluate  the  strains,  but  more 
easily  we  can  use  the  fact  that  the  4-node  shell  element  of  ref.  [3] 
satisfies  the  rigid  body  mode  criterion.  Hence,  for  a  rigid  body  dis¬ 
placement  the  anc*  ^23  are  zero»  ^ro[T1  which  it  follows  that  also 

the  shear  strains  in  Eq.  (3)  are  zero,  and  the  rigid  body  mode  criterion 
is  satisfied. 

2)  The  element  can  approximate  the  Kirchhoff-Love  hypothesis  of  negli¬ 
gible  shear  deformation  effects  and  can  be  used  for  thin  shells. 

Various  demonstrative  solutions  are  given  in  Section  4. 

3)  Based  on  our  studies  the  element  does  not  contain  any  spurious  zero 
energy  modes  (using  a  "full"  numerical  integration). 

We  reach  this  observation  by  studying  the  strains  along  the  element 
sides.  If  the  element  were  to  contain  a  spurious  zero  energy  mode,  the 
strains  along  every  side  should  vanish  for  a  displacement  pattern  (to  be 
identified)  other  than  the  displacements  corresponding  to  a  true  rigid 
body  mode.  However,  such  displacement  pattern  could  not  be  identified. 

Considering  the  practical  use  of  the  element  the  interpolation 
employed  for  the  transverse  shear  strains  shows  that  1S  constant 

with  r-j  and  in  general  discontinuous  at  r^  =  +_  1  (between  elements), 

and  similarly  is  constant  with  r^  and  in  general  discontinuous 

at  *2  =  —  1 •  As  a  consequence,  the  accuracy  with  which  transverse 

shear  stresses  are  predicted  depends  to  a  significant  degree  on  the 
mesh  used  and  the  geometric  distortions  of  the  elements.  However, 
our  experience  is  that  the  bending  stress  predictions  are  relatively 
little  affected  by  element  distortions  (see  Section  4). 

To  employ  Eq.  (7),  we  also  need  to  use  the  appropriate  constitutive 
relations: 


(8) 
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~i  iki 

where  C  is  the  fourth-order  contravariant  constitutive  tensor  in  the 
convected  coordinates  r..  The  constitutive  law  is  known  in  the  local 

Cartesian  system  of  orthonormal  base  vectors  e.,  i  =  1,  2,  3,  with  the 
*33  — 1 

condition  f  equal  to  zero,  see  Fig.  3  (refer  to  [2],  p.  258).  Denoting 

this  constitutive  tensor  by  Cmnop,  the  constitutive  tensor  for  Eq .  (8)  is 
obtained  using  the  transformation 


:i  jk;. 


(g1 


lm> 


(gv 


in)  <£' 


io>  (£ 


e p  )  Cmnop 


(9) 


where  the  are  the  contravariant  base  vectors  of  the  convected  coordi¬ 
nates  r. .  These  vectors  are  calculated  using  the  co variant  base  vectors 

£. ,  where 


(10) 


with  0^  from  Eq .  (1)  and  the  following  relations. 


and 


i 

g 


h 


on 


(12) 


where  D1^  is  the  cofactor  of  the  term  g. .  in  the  matrix  of  the  metric 

*  J 

tensor  and  |Jj  is  the  determinant  of  the  Jacobian  matrix  at  the  point 
considered. 
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?3 

-g3| 


*1 


h  x  -3 

S2  x  % 


*2 


-3  x  *1 


Figure  3  Local  Cartesian  coordinate  system  used 
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3.  TOTAL  LAGRANGIAN  FORMULATION 


The  large  displacement  formulation  of  the  shell  element  is  based  on 
the  derivation  given  in  ref.  [2,  Section  6.3.5],  and  the  concepts  and 
interpolations  presented  in  the  previous  section. 

The  geometry  of  the  element  at  any  time  t  is  defined  as  in  Eq.  (1) 


t  k 

but  using  the  nodal  point  coordinates,  x., 

(A)  1 
time  t,  ' 


t  k 

and  director  vectors  V  ,  at 


t 


(13) 


where  we  imply  summation  over  k.  The  displacements,  tu. ,  and  incremental 


displacements,  u. ,  of  a  particle  of  the  element  at  time  t  are  hence  given  by 


t 


(14) 


t  k  k 

where  the  u.  are  the  nodal  point  displacements  at  time  t,  the  u. 


are  the  incremental  nodal  point  displacements  from  the  configuration 


t  k  t  k 

at  time  t,  and  the  variables  ,  V^. ,  and  are  defined  as  in 


Eq.  (2)  but  referred  to  the  configuration  at  time  t. 

This  kinematic  description  implies  the  following  hypotheses: 

•The  director  vectors  remain  straight  during  the  deformations. 

•The  "thickness"  of  the  element  measured  along  the  director 
vectors  remains  constant  during  the  deformations;  hence  only 
small  strain  conditions  are  considered. 

Using  the  assumptions  in  Eqs.  (13)  and  (14)  the  geometric  and  material 
nonlinear  response  is  analyzed  using  an  incremental  formulation  [2],  in 
which  the  configuration  is  sought  for  time  (load  step)  "t  +  At",  when  the 
configuration  for  time  "t"  is  known.  The  basis  of  this  incremental 

^  ^Note  that  the  superscript  t  on  a  variable  denotes  the  configuration 
at  time  t  in  the  incremental  solution  and  does  not  imply  a  dynamic 
analysis . 
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formulation  is  the  use  of  the  virtual  work  principle  applied  to  the 
configuration  at  time  t+At.  In  essence,  two  approaches  can  be  employed 
leading  to  the  updated  Lagrangian  and  the  total  Lagrangian  formulations. 
These  approaches  are,  from  a  continuum  mechanics  point  of  view,  equivalent 
and  in  the  following  we  develop  the  governing  finite  element  relations  for 
the  total  Lagrangian  formulation. 

The  principle  of  virtual  work  applied  to  the  configuration  at  time 
t+At  is 


t+At  ~ i j 
05 


t+AtR 


(15) 


where  the  ^+“gS1J  are  the  contravariant  components  of  the  second 

Piol a-Kirchhoff  stress  tensor  at  time  t+At  and  referred  to  the  configurati 

at  time  0,  and  the  are  the  covariant  components  of  the  Green- 

Lagrange  strain  tensor  at  time  t+At  and  referred  to  time  0.  Both  sets  of 
tensor  components  are  measured  in  the  convected  coordinate  system  r . , 

4-1  i  f  I 

i  =  1,2,3.  The  external  virtual  work  is  given  by  atR  and  includes 
the  work  due  to  the  applied  surface  tractions  and  body  forces. 

For  the  incremental  solution,  the  stresses  and  strains  are  decomposed 
into  the  known  quantities,  gS1"^  and  and  unknown  increments,  gS1^ 

and  gE.j,  so  that 


t+At^ij 

=  + 

0 

0 

t+At~ 

0£ij 

U  , 

"  0£ij 

(16) 

(17) 


In  addition,  the  strain  increment  can  be  written  as  a  linear  part, 
ge^.,  and  a  nonlinear  part,  hence 
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(18) 


Substituting  from  Eqs.  (16)  to  (18)  into  Eq.  (15)  and  using  the 
linearized  expressions  0S1J  =  gC1^1  Qekl  and  =  50®ij»we  obtain 

the  linearized  equation  of  motion 


C*kl 


Ov 


+ 


a  t+AtR 


(19) 


This  equation  is  the  basic  equilibrium  relation  employed  to  develop  the 

governing  finite  element  matrices.  For  the  actual  solution  of  problems. 

it  is  frequently  important  to  use  equilibrium  iterations,  but  the  finite 

element  matrices  and  vectors  used  in  these  iterations  can  be  derived 

directly  from  the  matrices  obtained  using  Eq.  (19)  [2].  Note  that 

' ' w  1  i  1c  1  /  tz  33 

nC  J  is  now  obtained  using  Eq.  (9)  with  the  condition  nS  =  0,  which 

implies  the  more  natural  condition  t  =  0  only  in  the  small  strain  case. 

The  basic  problem  of  the  finite  element  discretization  of  Eq.  (19) 
lies  in  expressing  the  strain  terms  of  Eq.  (19)  in  terms  of  the  finite 
element  interpolations.  Using  the  definition  of  the  Green-Lagrange 
strain  components 


t. 

0Eij 


(20) 
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and  the  relations  in  Eqs.  (13)  and  (14)  we  obtain 

r-.  ^  a.  k 


0ei  i  -  hk 


V  Hk  +  Takhk  ,i  <-“k  'll  '  %  +  6k  **  -M)  <21-a) 


0ni  i 


y-  hk,i  hp,i  ik  •  ip  +  “rhk,1  hp,1  ap('“p  '  ik  +  8p  M  ■  Mfc) 


(r3^  .  .  ,  t..k  .  0  t.,kk 

+  T  hk,i  hp,i  ak  ap  ("ak  ^2  +  Bk  V 


(-a  +  B  )  (21-b) 


(i  =  1,2) 


9^k  T  k  k  k 

with  the  notation  hk  ^  ,  _uk  =  [u-|  u2  u-j],  and 


0e12 


2 


hk,2  lai  •  ak+  hk,i  •  ik 


u,  + 


r3  .  .  /  t„k  t_  j_  „  tuk  t. 


2  k  ,2  k 


/  t„K  L  ,  0  U,K  \  , 

(-ak  V2  .  +  8k  V,  .  )  + 


T  hk,1  ak  ('°k  ^—2  *  ^2  +.  e 


M  •  V] 


°"12  -  T[hk,l  hp,2  ik  “, 


(22-a) 


“f-  hk,l  hp,2  (-«p  %P  •  ufc  +  6p  Sf  •  uk)  + 

\,1  hp,2  ^(-“k^-ip^kM  •  V  + 

\,2  \  ap  K  **2  +  8k  *tf>  •  <-“p  *£  +  3P  ^1>. 


(22-b) 
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Further,  we  obtain  for  the  transverse  shear  strains,  using  Eqs.  (3) 
and  (6), 


oei3  =  ~r 


<1  *  r;)  [  V,  („|  -  a,’^  •  6,  a,  Vu 

M,i] 


V2  tv2i  +  62a2 


+  T  (1  "  r2)  [tg3i  (ui  '  ui3)  +  T  tglCi  ("a4a4  Mi  +  64  a4 


'a3  a3  ty2i  +  B3  a3 


(23-a) 


°r13  =  32  0  +  r2)  [('al  ai  t'-*2i  +  B1  al  tvii  “  a2  a2  ty2i  +  32  a2  Mi^ 

(u!  - 

+  32  d'  r2^fr"a4a4  V2i  +  B4a4  ^li  '  a3  a3  tv2i  +  B3  a3tyi i ^  (ui  '  ui^] 


(23-b) 


and, 


0  23  8 


r  0+n)  [l931  (u!  -  u?)  +T  921  (-Vl  V21  +  6lal  *v]i 


°4a4  ty2i  *  B4  a4  Ml>] 


T  (1-r,)  [‘g®,  (u*  -  u?)  %  (-V2  Ml  +  «2a2  tv 


t„2 
li 


’  a3a3  tv2i  +  B3a3  Mii 


(24-a) 


21 


nt 


i 


U  23  32 


(1  +  r-] ) 


Vl  tV2i  +  8lal  tv!i  *  “4a4  tv2i  +  84a4  tvii>  <ui  -  u1>] 


32  H  -rj)  ^(-o2a2  V21  +  fi2a2  - 


a3a3  ‘V2i  +  83a3  tvii)  <"l 


(24-b) 


-  «?>] 


Note  that,  since  we  assume  the  thickness  of  the  shell  to  be  constant, 
the  strain  through  the  element  thickness  is  zero. 

The  expressions  in  Eqs .  (21)  to  (24)  are  substituted  into  Eq.  (19) 
which  in  the  standard  manner  yields  the  linear  strain  incremental  stiffness 

matrix  the  nonlinear  strain  (or  geometric)  incremental  stiffness  matrix 

t  t 

and  the  nodal  point  force  vector  in  the  finite  element  incremental 

equilibrium  relations  [2], 


( +  o%  ^ 


u  = 


t+Atp  t, 

-  o- 


(25) 


The  element  matrices  in  Eq.  (25)  correspond  to  five  degrees  of 
freedom  Der  node,  see  Fig.  1,  but  in  some  applications  it  is  convenient 
to  use  instead  of  and  8^  three  rotations  about  the  global  coordinate 

axes  (see  Section  4.6).  In  this  case,  we  simply  transform  the  matrices 
of  Eq.  (25)  in  the  standard  manner  [2]. 

The  element  developed  can  be  used  for  nonlinear  dynamic  analysis  of 
shells  with  any  time  integration  method  [2],  using  a  consistent  mass  matrix 
calculated  in  the  standard  way  or  a  lumped  mass  matrix  with  the  diagonal  terms 

°v  0  °v  2 

m^  =  — c-^  -  and  m. .  =  (a^)  f°r  translational  and  rotational  degrees 

of  freedom,  respectively. 
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4.  NUMERICAL  TESTS  AND  EXAMPLE  SOLUTIONS 


We  have  implemented  our  shell  element  in  the  ADINA  computer  program 
and  have  performed  various  numerical  tests  to  study  the  predictive  capa¬ 
bilities  of  the  element.  The  following  solutions  were  all  obtained  using 
2x2  Gauss  integration  in  the  r3  =  0  surface  of  the  element,  and  2  and  4 

point  Gauss  integration  in  the  r^  direction,  for  elastic  and  elasto-plastic 

analyses,  respectively. 

4.1  Some  simple  tests 

As  a  first  step  to  test  the  element,  the  eigenvalues  of  the  stiffness 
matrices  of  undistorted  and  distorted  elements  were  calculated.  In  all  cases, 
as  expected,  the  element  displayed  the  six  rigid  body  modes  and  no  spurious 
zero  energy  modes. 

4. 1 .a  Patch  tests 


For  the  patch  test  [2,  18]  the  mesh  shown  in  Fig.  4  was  used.  In  the 
first  analysis  the  mesh  was  loaded  with  the  constant  moment  indicated,  and  a 
constant  curvature  (linear  distribution  of  rotations)  was  obtained  for  both 
plate  thickness  in  the  two  plate  directions.  The  transverse  displacements 
predicted  by  the  model  were  —  as  expected  —  those  of  Ki rchhoff-Love  plate 
theory  at  nodes  7  and  8. 

In  the  second  analysis  the  rotational  degrees  of  freedom  were  deleted 
and  the  mesh  was  subjected  to  shear  forces.  As  expected,  for  both  plate 
thicknesses  a  linear  distribution  of  transverse  displacements  was  obtained. 

In  the  third  analysis  the  mesh  was  subjected  to  an  external  twisting 
moment.  In  the  thin  plate  analysis,  constant  curvatures  were  obtained  in 
both  plate  di  rections,  and  the  transverse  displacements  agreed  with  the 
analytical  thin  plate  solution.  In  the  thick  plate  analysis,  a  slight 
non-symmetry  in  the  displacement  response  (the  third  digit)  was  obtained 
due  to  the  unsymmetric  representation  of  the  transverse  shear  deformations. 
This  non-symmetry  is  not  observed  if  the  shear  deformations  are  suppressed 
(which  corresponds  to  thin  plate  theory)  by  choosing  a  large  value  for  the 
shear  correction  factor  k  [2,  p.  236]  (or  when  using  rectangular  elements 
i n  the  mesh) . 

Finally,  it  should  be  noted  that  the  patch  test  is  of  course  passed 
for  the  three  membrane  stress  states  (t^,  and  t-^  constants). 

4 . 1 . b  Cantilever  linear  analyses 


A  cantilever  of  unit  width,  thickness  0.1  and  lengths  10  and  100 
was  subjected  to  a  tip  bending  moment.  The  structure  was  modeled  using  one 
single  element  and  two  distorted  elements  as  shown  in  Fig.  5.  The  results 
obtained  in  these  analyses  for  the  displacements  and  rotations  ct  the 
cantilever  tip  and  the  stresses  were  those  of  Bernoulli  beam  theory. 
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E  =  2.1  x  106 


a)  One  element  case 


Figure  5  Cantilever  subjected  to  tip  bending  moment 
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Next,  the  cantilever  in  Fig.  6(a)  was  analyzed  for  the  transverse 
tip  load  shown.  Using  4  equal  size  elements  to  idealize  the  cantilever, 
again  good  results  were  obtained  when  compared  with  beam  theoretical 
resul  ts . 

Finally,  the  elements  modeling  the  cantilever  were  distorted  as  shown 
in  Fig.  6(b)  for  a  thin  and  a  thick  cantilever.  The  results  show  that 
the  transverse  displacements  and  normal  bending  stresses  are  almost 
insensitive  to  the  element  distortions.  However,  the  calculated  trans¬ 
verse  shear  stresses  (not  shown  in  the  figure)  are  not  accurate. 

4.1.c  Linear  analyses  of  a  simply-supported  plate 

A  simply-supported  plate  was  considered  for  a  static  and  a  fre¬ 
quency  analysis  using  a  consistent  mass  matrix.  To  model  one  quarter 
of  the  plate  a  4x4  mesh  of  equal  elements  was  used.  Figure  7  gives 
a  comparison  of  the  numerically  and  analytically  predicted  results.  The 
same  plate  was  also  analyzed  using  the  distorted  element  mesh  also  shown 
in  Fig.  7(a)  and  the  results  of  Fig.  7(c)  were  obtained. 

4 . 1  . d  Analysis  of  a  rhombic  cantilever 

The  rhombic  cantilever  shown  in  Fig.  8,  fixed  at  one  side  and  sub¬ 
jected  to  constant  pressure  was  analyzed  using  a  4x4  element  mesh.  The 
results  for  the  transverse  displacements  at  six  locations  are  compared 
against  the  solutions  obtained  using  the  DKT  triangular  element  of 
ref.  [6],  experimental  measurements  [1]  and  using  the  16-node  isopara¬ 
metric  element  (with  4x4x2  Gauss  integration).  In  all  cases  a  one 
step  geometric  nonlinear  analysis  with  equilibrium  iterations  was  per¬ 
formed.  Good  correspondence  between  the  experimental  results  and  the 
solution  obtained  using  our  new  4-node  element  is  observed. 

4 . 2  Linear  analysis  of  a  cylindrical  (Scordeli s-Lo)  shell 

The  shell  structure  shown  in  Fig.  9  has  frequently  been  used  to 
test  the  performance  of  shell  elements  [12].  Figure  9  shows  the  solutions 
obtained  with  our  elements.  In  each  of  the  solutions  uniform  meshes  with 
equal  sized  elements  were  employed  over  one  quarter  of  the  shell.  Solutions 
obtained  using  the  3-node  DKT  triangular  element  [25]  and  the  16-node  iso¬ 
parametric  element  [25]  are  also  shown. 

4.3  Linear  analysis  of  a  pinched  cylinder 

The  pinched  cylinder  problem  shown  in  Fig.  10  was  also  frequently 
analyzed  to  test  shell  elements.  Figure  10  shows  the  convergence  behaviour 
obtained  with  our  new  element,  when  comparing  the  finite  element  solutions 
with  the  solution  given  in  refs.  [11,  21].  Note  that  using  the  isoparametri 
shell  element  of  ref.  [3],  also  a  fairly  large  number  of  degrees  of  freedom 
are  required  to  predict  the  response  of  the  cylinder  accurately. 
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a)  Solution  using  non-distorted  elements 

Figure  6  Response  of  a  cantilever  subjected  to  transverse  tip  load, 
stresses  shown  are  those  at  the  Gauss  integration  stations 
r3  =  0.57735;  Tpp,  is  the  principal  stress  in  the  distorted 
mesh,  and  its  direction  was  always  less  than  11  degrees  from 

the  X2  axis. 
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b)  Solution  using  distorted  elements 


Figure  6  continued 
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Mode  shape 
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1-1 
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1-3 

1.18 

3-3 

1  .17 

b)  Natural  frequencies  (cycles/sec)  calculated  using  the  non-distorted  mesh 


Figure  7  Linear  analysis  ot  a  simply-supported  square  plate,  the  parameter 
of  distortion,  a,  was  equal  to  2.50. 
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c)  Static  response  due  to  constant  pressure  loading,  stresses 
are  given  along  line  Xg  =  0,  x^  =  0.028868. 

Figure  7  continued 
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Figure  8  Response  of  rhombic  cantilever  subjected  to  constant  pressure 
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R  =  300. 

L  =  600. 

*  =  40° 


Figure  9  Linear  analysis  of  a  cylindrical  shell  subjected  to  dead  weight. 

The  2x1  result  refers  to  the  solution  obtained  with  two  16-node 
shell  elements  spanning  from  C  to  B.  The  16x16  result  refers  to 
the  use  of  512  equal  triangular  DKT  elements. 
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b)  Comparison  between  4-node  and  16-node  elements 

Figure  10  Linear  analysis  of  a  pinched  cylinder;  u  =  axial 
displacement,  w  =  radial  displacement 
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c)  Displacements 
Figure  10  continued 
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4.4  Large  deflection  analysis  of  a  cantilever 


The  cantilever  shown  in  Fig.  11  was  analyzed  for  its  large  displace¬ 
ment  and  large  rotation  response.  This  is  a  typical  problem  considered 
to  test  the  geometric  nonlinear  behaviour  of  beam  and  shell  elements  [25]. 
Figure  11  shows  also  the  models  used  in  the  analysis. 

The  first  two  models  are  single  element,  cubic  and  parabolic  iso¬ 
parametric  degenerate  shell  element  models.  Model  I  predicts  the  response 
of  the  cantilever  very  accurately,  whereas  model  II  yields  an  accurate 
response  solution  in  linear  analysis  but  locks  once  the  element  is  curved 
in  the  nonlinear  response  solution.  This  observation  is  in  accordance  with 
the  results  reported  in  ref.  [5]. 

The  same  nodal  point  layouts  were  next  employed  for  models  III  and 
IV  using  our  new  4-node  shell  element.  Figure  11  gives  also  the  results 
obtained  with  these  models.  It  is  seen  that  model  III  yields  an  accurate 
large  displacement  response  prediction,  and  even  model  IV  yields  quite 
accurate  results  up  to  about  60  degrees  of  rotation.  The  computer  time 
required  in  these  analyses  were  only  different  using  models  I,  III  and 
IV. 


Another  important  result  is  shown  in  Fig.  12.  As  reported  in  ref.  [5], 
the  cubic  shell  element  is  sensitive  to  "in-plane"  distortions;  and  hence, 
it  is  interesting  to  study  the  effect  of  using  a  distorted  element  mesh 
in  the  analysis  of  the  cantilever.  Figure  12  summarizes  the  results 
obtained  using  the  one  cubic  element  and  three  4-node  elements  with  a 
nodal  layout  that  corresponds  to  distorting  the  elements.  It  is  seen  that 
the  predictive  capability  of  our  new  4-node  element  is  considerably  less 
sensitive  to  the  element  distortions. 

4.5  Geometric  nonlinear  response  of  a  shallow  spherical  shell 


Figure  13  shows  the  spherical  shell  that  was  also  analyzed  in  ref.  [3] 
with  one  cubic  shell  element,  modeling  one  quarter  of  the  shell.  To  test 
our  new  4-node  shell  element,  the  same  nodal  point  layout  as  in  ref.  [3] 
was  used,  giving  a  mesh  of  nine  elements.  Figure  13  shows  the  response 
calculated,  including  the  post-buckling  response  (not  reported  in  ref.  [3]) 
with  the  automatic  load  steDDing  algorithm  of  ref.  [4].  Good  correspondence 
with  the  analytical  solution  of  Leicester  [20]  and  the  solution  of  Horrigmoe 
[16]  was  obtained.  The  solution  with  the  16-node  element  was  almost  twice 
as  expensive  as  the  4-node  element  solution  (using  in  both  cases  the  same 
parameters  for  the  automatic  step-by-step  solution  algorithm). 

4.6  Linear  buckling  analysis  and  large  deflection  response  of  a  simply- 

supported  stiffened  plate 

The  stiffened  d1  ate  shown  in  Fig.  14  was  analyzed  for  its  buckling 
response.  Since  we  expect  the  buckling  mode  to  be  symmetric  [26]  only  one 
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a)  Finite  element  models 

Figure  11  Large  deflection  analysis  of  a  cantilever 
using  non-distorted  elements 
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b)  Response  of  model  I 
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c)  Response  of  model  III 
Figure  11  continued 
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d)  Response  of  model  IV 


Figure  11  continued 
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Figure  12  Large  deflection  analysis  of  a  cantilever  usinq  distorted 
elements 
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Figure  13  Geometric  nonlinear  response  of  a  spherical  shell 
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<L 


simply  supported  plate 
E  =  2.1  x  106 
v.  =  0.3 


Figure  14  Nonlinear  response  of  a  stiffened  plate 
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quarter  of  the  plate  is  modeled  using  symmetry  boundary  conditions.  The 
model  consists  of  nine  4-node  shell  elements  and  three  2-node  isoparametric 
beam  elements.  At  the  nodes  where  a  shell  element  connects  to  a  beam 
element,  three  rotational  degrees  of  freedom  aligned  with  the  global  axes 
are  considered  for  the  shell  element.  In  order  to  avoid  locking  of  the 
isoparametric  beam  elements,  one  point  Gauss  integration  along  the  beam 
axes  was  used.  This  does  not  introduce  spurious  zero  energy  modes  in  the 
model,  although  the  bending  stiffness  of  the  beam  is  under-estimated. 

The  linearized  buckling  problem  was  solved  as  described  in  [4,  Eq.  (37)] 
and  we  obtained 


°cr  element  solution) 

0cr  (analytical  solution) 


Next,  an  initial  imperfection  with  the  shape  of  the  first  buckling 
mode  and  a  maximum  amplitude  of  1 / 5th  of  the  plate  thickness  was  intro¬ 
duced.  Figure  14  shows  the  large  deflection  response  of  this  model  as 
calculated  using  the  automatic  load  stepping  scheme  of  ref.  [4]  with  a 
tight  enerqy  convergence  tolerance. 

4.7  Analysis  of  el asto-olastic  response  of  a  circular  plate 

The  thin  circular  plate  shown  in  Fig.  15  was  analyzed  for  its 
elasto-plastic  response,  when  subjected  to  a  concentrated  load  at  its 
center.  The  plate  is  simply-supported  with  its  edges  restrained  from 
moving  in  its  plane. 

In  a  first  solution,  the  plate  model  shown  in  Fig.  15(a)  was 
used  to  analyze  the  plate  assuming  small  displacements  (materially- 
nonl i near-only  conditions).  Figure  15  shows  that  the  theoretical 
collapse  load  is  over-estimated,  but  for  the  coarse  mesh  used,  the 
predicted  response  is  quite  reasonable. 

In  a  second  solution,  large  displacements  and  elasto-plastic 
conditions  were  assumed  and  in  this  case  the  stiffening  behaviour  of 
the  plate  shown  in  Fig.  15(b)  was  predicted.  In  order  to  have  a  com¬ 
parison,  also  the  model  of  five  axis.ymmetric  8-node  elements  shown  in 
Fig.  15(a)  was  solved.  Fioure  15  shows  that  both  models  predict  in 
essence  the  same  response;  however,  in  this  case  relatively  little 
plasticity  was  developed  for  the  range  of  displacements  considered. 

4.8  Dynamic  analysis  of  an  el asti c-perfectly  plastic,  simply-supported 

square  plate 

One  quarter  of  the  plate  shown  in  Fig.  16  was  modeled  using  sixteen 
4-node  elements.  The  central  difference  method  was  used  in  the  time 
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R 


ax i symmetric  model 
a)  Finite  element  models 

Figure  15  Response  of  elastic-perfectly  plastic  circular  plate 
subjected  to  a  concentrated  load,  P,  at  its  center. 
TLF  abbreviates  use  of  total  Lagrangian  formulation 
and  MNO  abbreviates  use  of  material ly-nonl inear-only 
formulation. 


106  et  =  0.0 
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b)  Circular  plate  response 


Figure  15  continued 
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b)  Time  hi  story 

Fio.  16  -  Dynamic  Analysis  of  Elastic-Plastic  Plate 
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integration  to  calculate  the  time  history  response  when  the  plate  was 
subjected  to  a  step  uniform  pressure.  A  lumped  mass  discretization 
was  employed. 

In  order  to  assess  the  effect  of  the  lumped  mass  matrix  on  the 
integration  step  (At),  the  smallest  period  of  the  finite  element  model 
was  calculated  using  a  consistent  mass  matrix  and  a  lumped  mass  matrix; 
we  obtained 


JUMPED 

jnin _  =11. 

tconsist. 

min 

For  the  transient  solution  we  used  the  material-nonlinear-only 

formulation  and  At  =  2x10"^.  Figure  16  compares  oredicted  vertical 
displacement  of  the  center  of  the  plate  with  the  solutions  given  in 
[3]. 


5.  CONCLUSIONS 


A  new  four-node  non-flat  general  nonlinear  shell  element  has  been 
presented  with  the  following  important  element  properties: 

-  The  element  is  formulated  using  three-dimensional  continuum 
mechanics  theory;  hence,  the  use  of  the  element  is  not  restricted 
by  application  of  a  specific  shell  theory. 

-  The  element  is  reliable  and  has  good  predictive  capability  in 
the  analysis  of  thick  and  thin  shells  . 

-  The  amount  of  computations  required  to  calculate  the  element 
stiffness  matrix  are  very  closely  those  that  are  used  in  standard 
isoparametric  formulations.  The  computer  time  used  could  be 
reduced  considerably  in  elastic  analysis  by  using  analytical 
integration  through  the  element  thickness. 

In  this  paper  we  have  presented  the  formulation  and  some  appli¬ 
cations  of  the  element.  The  solution  results  obtained  are  most  encouraging, 
but  a  formal  mathematical  convergence  study  of  the  element  would  be  very 
valuable,  and  we  are  currently  pursuing  such  research. 

Finally,  it  should  be  noted  that  the  element  presented  here  provides 
a  very  attractive  basic  formulation  that  could  be  extended  to  large  strain 
analysis  and  analysis  of  composite  shells.  Also,  the  concepts  applied  here 
to  formulate  a  4-node  element  could  equally  well  be  employed  in  an  effective 
manner  to  formulate  higher-order  shell  elements. 
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